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INTRODUCTION 

Similarities between the statistical behavior of molecules in a gas and 

the velocity fluctuations of fluid elements in a turbulent flow suggest the 

possibility of describing both phenomena in terms of a velocity distribution 

function from which mean properties may be computed by taking appropriate 

moments. The literature abounds with efforts to exploit this analogy between 

the two areas of the field of statistical mechanics. From another point of view, 

there are potential simplifications in modeling turbulence behavior at the 

level of the velocity distribution function (or, the probability density) , 

rather than modeling Individual correlations. Of particiilar interest are two 

recent studies which have attempted to develop theoretical treatments for 

12 1 

practical application ’ . Lundgren begins with the incompressible Navier- 
Stokes equations in a form which assumes either an infinite fluid region or a 
constant pressure boundary condition. A hierarchy of equations for multipoint 
distribution functions is developed which strongly resembles the Bogoliubov- 
Bom-Green-Kirkwood-Yvon^ (BBGKY) equations. In a subsequent work Lundgren^ 
attempts to close the system at the one-point level by employing a relaxation 
model identical in form to the Bhatnager-Gross-Kfook (BGK) model of kinetic 
theory. This model equation is not, within Itself, sufficient to define a 
turbulent flow. An additional equation is required to relate the turbulence 
dissipation rate to other flow properties. This, in effect, implies that an 
ad hoc assumption must be made regarding the relaxation rate in the model. 
Lundgren applies his eqmtions to several idealized probl ems in which no solid 

boundaries are present. 

2 6 

Chung ’ has taken a very different,, approach but arrived at a similar 
governing equation for the distribution function. His analysis was developed 
from ideas of generalized Brownian motion and resulted in a modified 


Fokker-Planck equation. The analysis has been extended to account for chemical 

7 

species and reactions . Chung likewise finds it necessary to malce M hoc 
assumptions on mixing length if his equations are to be self-contained. The 
resulting model is applied to the problem of plane Couette flow of an incom- 
pressible, single-species gas by employing moment methods familiar in kinetic 

g 

theory . In these methods specific fimctional forms are assumed for the dis- 
tribution function and unknown coefficients are found in the solution. 

The present work describes the solution of Lundgren’s model equation 
for plane Couette flow. This provides an important extension of his previous 
studies in that the flow field is boimded by solid surfaces and in that it 

represents a flow for which experimental data are available. The solution 

9 10 

is acconplished by an extension of the discrete ordinate method ’ as 
developed for problems in rarefied gasdynamics. This differs from the moment 
method selected by Chung in that no a priori assumption is enforced upon the 
form of the distribution function. Res\alts obtained in the present study are 
compared with Chung’s work and with available experimental data. 

Perhaps an interesting analogy between the application of the discrete 
ordinate method and the moment method is that of the improvement provided in 
turbulent boundary layer calculations by the use of numerical solution to the 
partial differential equations themselves as opposed to the classical integral 
method used earlier. A finer detail of flow structure is afforded by the 
direct numerical solution. 


2 



GOVEEDUKG EQUATIOMS 


General Equation for the Distribution Function 

The starting point for the present analysis is the lowest order equation 
for the turbulent distribution function (Eqs. (l) and (2), Ref. k) 


■3f 

3t 



f 


dv 


(1) 


where f(r, v, t) dv is the probability that the velocity at point r in 
physical space is in the range v to v + dv in velocity space. Pressure, 
density and kinematic viscosity are symbolized by p, P and v , respectively. 
The relaxation time is denoted by T, which is related to the characteristic 
turbulence diffusion time. F is a Gaussian (equilibrium) distribution given 
ty 


= ^2n' 




[- (v - u)^2U^ 


( 2 ) 


The time average flow velocity is u and 3bT is the mean square of the 
velocity fluctuations. These are defined as appropriate moments of f. That is, 
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= J 


3 

fd^v 


3Ij2 = J (tt _ u)‘ 
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fd-^v 


where the integrations are taken over the entire velocity space ( - “ to + “ 
for each component) . 
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Ltuidgren assumes that the relaxation time is approximately l/u , where L 
is the integral scale of turbulence, and models this by the equation 


K(e + I DU^/Dt) 


(3) 


Here, e is the turbulence dissipation rate, D/Dt is the usual substantial 
derivative, and K is taken to be a constant whose value is approximately 5- 
If the turbulence velocity c = v - u is introduced as an independent 
variable, the analysis is made somewhat simpler. Making this change of 
variable and simplifying the physical coordinates to the case of one dimension 
in anticipation of the application to Couette flow, one obtains 
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Here, y is taken to be normal to the mean flow and 
turbulence velocity components. 
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Reduced Distribution Functions 

A reduction in computer storage requirements is afforded by defining re- 
duced distribution functions according to the following scheme: 


Cy) = J V ""y’ ‘^''x 
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Also, let 
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(^’ ^y) ^ I I (<^X ■" 0 =y’ ^^z) ^'^x 

••CO 

The reduced equilibrium distribution fxmctions, G, J, H and J are defined 
as shown above with f = S’. Equation (4) may now be transformed into a system 
of equations which, although similar in form to the equation for f , are much 
easier to treat numerically. These become 
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These distribution functions should, satisfy the constraints 

CO 

J g dCy = 1 (7a) 

..CO 

and. “ 

J d = 0 (7b) 

..CO 

The first of these states that the probability of finding a fluid element 
somewhere in r , v space is unity, while the second requires that the mean of 
the longitudinal fluctuating velocity component should be zero. 

Moments of Interest 

Once the equations are solved for the distribution functions, moments 
of interest may then be calculated. For example, one obtains 


u = j dc 
X d "'v y 
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These are the mean -velocity, t-urhulence kinetic energy, Eeynolds stress, 
mean square of y-velocity component fluctuation, and kinetic energy flux, 
respectively. 

Final Reduced Equations 

It is convenient to define the following nondimensional variables : 
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where 2d is the distance between the plates, p is the turbulence shear _i 

Xy r -|2 

stress, is the kinematic viscosity, p is the density, and u^ = L^^xy^w' 
is the usual friction velocity. 

Using these quantities and dropping the superscript " " , for 

simplicity, the set of governing equations becomes 
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It is this set of equations which is to he solved by the discrete ordinate method, 
subject to the boundary conditions to be discussed subsequently. 
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Examination of this system of equations reveals that it is not yet 
completely self-contained. The dissipation rate e is not given as a moment 
of the distribution function. Therefore, a separate equation for e is 
required to close the set. Among the possibilities for such an equation for 
the present problem are the use of an algebraic relationship resulting from 
setting turbTolence production equal to dissipation 


e 



or 


= 0.7968 U^/y 


(9e) 


or by use of a differential equation to model e , such as the one derived by 
Jones and Launder^ based upon a semi-en 5 )irical approach. For Couette flow it is 
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’ e ^ 1 2 

Each of these equations for 

results will be discussed in 


constants and v = Ee„ 
1 * 

e has been employed in 
a subsequent chapter. 
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the present study, and 
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BOUHDAEY COMDITIOK8 


One of the difficult aspects of this research is specification of the 
appropriate boundary conditions. This difficulty arises as a consequence of two 
factors. First, since the mean velocity is a moment of the distribution function, 
information from continuum flow only gives the no-sllp condition that 


u 


i 


V f 


d\ = 


Because many functions for f could be specified which satisfy this integral 
constraint, there is a lack of uniqueness in prescribing f . If f is prescribed 
as a Dirac delta function which also in5)lles zero instantaneous velocity at the 
wall, it is difficult to Incorporate numerically. 

Second, the relaxation time, t , is difficult to model in the viscous sub- 
layer because of the unknown variation of e and U. Furthermore, in the 
derivation of the turbulence model equation, the effects due to the presence of 
the wall were not modeled separately. As a consequence, the validity of the 
model equation is suspected very near the wall even if a successful relaxation 
model is obtained. Thus, there is a hesitation in applying the analysis in this 
zone . 

With these considerations in mind there are clearly two fundamental 
questions to be addressed in the research. 

(i) where should the boundary condition on the distribution 
function be applied, and 

(ii) what functional form should it take? 

Thus, the present study can actually be divided into two segments: first, 

to develop the capability of obtaining convergent, stable numerical solutions to 

the model equation; and second, to examine the validity of various boundary 

12 1 ^ 

conditions in light of conpiarlson with experimental data ’ . This division of 

the problem virtually parallels the situation arising in calculations of rarefied 
flows from the Boltzmann equation. In this field the study of gas-STzrface 
interactions-- which are used for bo-undary conditions in solution of the govern- 
ing equation — has become almost a separate field of study within itself. Such 
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could well tie the case with the statistical approach to turbulent flows. 

Matching to law of Wall 

The course taken in the present work was to co3xfine the application of 
the turbulent model equation to regions outside the viscous sublayer. Values 
of y* = from. 50 to 100 were selected as the boundary points for the 

governing equation, and the usual functional forms for law of the wall were 
assumed to relate the boundary point to wall conditions. This necessitates 
that certain matching of the numerical solution to law of wall variation at 
the boundary point be performed. The specific conditions applied are dependent 
upon the form selected for the distribution function at the boundary and will 
be discussed in more detail in the section on MJMERICA.L APPROACH. 

Zero Gradient Distribution Function 

The momentum equation for Couette flow with zero pressure gradient in the 
continuum theory for turbulence reduces to a statement that total stress is 
constant between the plates. If attention is confined to the region well out- 
side the viscous sublayer, then the viscous stress contribution is negligible 
and therefore Reynolds stress is constant. Assuming that the apparent viscosity 
coefficient is linear in y then gives the familiar logarithmic result for the 

mean velocity profile. Further, the kinetic energy of turbulence, , is 

12 

known to be approximately constant in such a logarithmic region . In terms of 
the moments described earlier, these conditions give 

CO 

rrr fc^ + c^ + c^) f dc dc dc = 3U^ = constant 
JJJ \ X y z/ X y z 


CO 

JJJ ‘'x °y ^ ^‘^y = constant 

mmCO 


Therefore, under the assumption of a linear variation in for the 

Couette flow problem, it should be possible to construct a boundary condition 
for the distribution fixnctlon which causes the governing statistical equation 
to yield a logarithmic solution between the plates. A necessary condition for 
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this is to forbid f to have a gradient at the boundary point and to require 
it to match the law of the wall there. 

The governing equation itself may be used to develop an appropriate form 
for f by setting 9f/9y = 0 in Eq. (4). If one then assumes du/dy = l/>ty, 
U - U_, e = du/dy , and 1/t = K e/u^ , the equation becomes 

i3 


3 =y ^ ar = 3 K(^b - " 3fg + ^ 

X X 


( 10 ) 


G 5 + C T 

y oc z oc 

y z 


where the pressure gradient and viscous stress have been taken as zero-. The 
subscript B is used to indicate that the distribution function obtained from 
solving this equation is to be employed as a boundary condition for solution 
to Eq. (4) for the region between the plates. 

Again, it is found to be easier to work with the reduced distribution 
functions. The equations for these which correspond to Eq. (lO) are 


3K(Gb 
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Solutions to these equations are independent of y 
"boundary condition functions in solving Egs. (6). 

0 since the mean velocity is not constant 
3 can he related to g and j hy 


and are to he applied as 
It should he noted that 
, hut logarithmic. However 
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n c f dc dc 
X X z 


+ u 


n f dc dc 
X z 


or 


= J + u g 


( 12 ) 


Hence, at the boundary 



% 


where Ug 


is determined hy the logarithmic relation 


u^ 





+ Bu^ 


with H and B as constants. A value for y^ = u^y/v of 100 is used to 
insure that viscous stresses are, in fact, negligible at the boundary point. 

The solution to the set of equations (ll) must he obtained numerically. 
Since the distribution functions and all their velocity gradients must approach 
zero as | c^l ^ this is used as a boundary condition an velocity space. 

The integration proceeds from a large value of c^ toward zero and from a small 
negative value c (with large absolute value) toward zero. A first order 
finite difference scheme is utilized for the integration. The result of the 
integration is a set of numerical values for fg which is then en^jloyed as 
a houndaay condition an the solution of the model equation for y^ > 100. 
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Chapman-Enskog Distribution Function 

There are two shortcomings in the use of the 9f/Sy = 0 boundary condi- 
tion. One is that the value of u^ must he specified a priori . This is 
equivalent to establishing the wall shear stress, which is a quantity which 
hopefully would emerge from the solution rather than be required as an input 
in order to obtain a solution. The second is that experimental data for mean 
velocity may not follow the logarithmic variation throughout the entire region 
between solution boundaries. This point will be discussed later in the 
EESULTS section. 

Therefore, it was desirable to seek a boundary condition for the distribu- 
tion function which would circumvent these difficulties. Upon first inspection 
it would seem possible to impose a Gaussian distribution for f^ in a manner 
analogous to the use of Maxwellian re -emission of molecules from a surface in 
the kinetic theory of gases. However, since a Gaussian distribution function 
gives zero Reynolds stress, this is inappropriate for application within a 
turbulent zone. (There may be merit in attempting to apply the statistical 
model equation within the viscous sublayer, where Reynolds stresses are small 
and then imposing the Gaussian distribution in a limiting Dirac delta function 
form at some point in this region. However, this possibility was not examined 
under the present effort) . 

An organized manner of obtaining a proper bounda:ry condition is the 

l4 

Chapman-Enskog procedure . Employing this method one can obtain approximate 
solutions to Eq. (l) using a series expansion. The zeroth order solution gives 
an equilibrium Gaussian distribution, which, as mentioned previously, results 
in zero Reynolds stress. The first order solution is commonly termed the 
Chapman-Enskog distribution, and allows for a Reynolds stress to occur. 

The first order Chapman-Enskog expansion for a one dimensional flow gives 


= p 


- f 
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dU^ + to “j\ 




dy dy JJ 


(13) 


where 
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2^2,2 
e + c + c 




2U^ 


is the Gaussian distribution. The corresponding nondimensional forms for the 
reduced distribution functions are 


gd) = <3 . ^ ffl c 
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These forms may then be applied as boundary conditions for the governing 
equation. Details of the application will be discussed in the chapter on 
RUMEEICA.L APPROACH. 

Two-Stream nature of Boundary Conditions 

Even though a functional form is established for the boundary distribution 
function, f^, , the correct implementation of this form is not straightforward. 
If one examines the physics of the problem, it is clear that both plates contri- 
bute to the establishment of the flow; and, therefore, boundary conditions 
should be applied at a value near each plate, giving two boundary condi- 

tions. (The symmetry condition at the centerline may be invoked to reduce 
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the problem to a half-space in y, but this still necessitates specifying two 
boundary conditions on f ) . Yet, if one examines the one -dimensional governing 
equation (Eq. (l)), it is observed that only a first-order derivative in y 
appears when viscous terms are neglected. It would thus appear that imposition 
of two boundary conditions would result in overspecifying the problem. 

Experience in solving the Boltzmann equation in rarefied gas dy^s^iics 
gives insight to resolving this paradox. In the molecular approach to rarefied 
flows one can only specify the velocity distribution function of molecules 
leaving a surface. The distribution function for those striking the surface is 
determined as a consequence of the solution. Thus, the boundary conditions 
possess a "two-stream" nature. The interaction of the incoming and outgoing 
stream is controlled through the collision or relaxation term in the model 
equation and through integral constraints such as the requirement that the 
incoming mass flux equal that for the outgoing stream. 

If this concept is transposed to the present problem of turbulent Couette 
flow, one requires that at the boundary point near the lower plate the distri- 
bution function is specified only for positive values of c while for the 

y 

corresponding point near the upper plate it is specified only for negative 
values of c . This is illustrated in Figure 1. Insofar as the function f 

y 

is concerned, this is equivalent to imposing a single constraint for all c^ 
domain values while it allows the effects of each plate to be introduced into the 
problem. This is mathematically consistent with the first order nature of the 
y-derlvative in the governing equation. Further, it seems plausible that such 
a two-stream approach is justified on a physical basis, since the turbulence 
motions leaving and approaching the wall region will be affected differently 
by the presence of the wall. 

As with the case in rarefied gasdynamics. Integral constraints must be 
imposed. The most obvious one is that 


CO 

f dc dc dc = 1 
JJJ X y z 

»co 


since this must hold from the definition of f. Additional constraints will be 
discussed in the chapter on NUMERICAL APPROACH. 
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3SIUMERICAL APPROACH 

Discrete Ordinate Method 

The discrete ordinate method is a ntmerical technique of replacing a 

continuous independent variable in a system of integro-partial differential 

equations by a set of discrete values and then treating these as parameters in 

the remaining solution. Although not restricted to dntegro-differentlal 

equations, the method has proven quite useful in attacking this type of problem. 

15 

Two esoamples of this application in physics are radiative heat transfer and 

9~10 

rarefied gas dynamics . The latter field is closely related to the present 
study of turbulence since a fundamental equation in rarefied gas dynamics is 
the Boltzmann equation for the velocity distribution function of molecules. 

If the BGK model is substituted for the collision integral of that equation, 
the one -dimensional form (which would correspond to a Couette flow geometry, for 
example) becomes 


c 1^ = - 7 (P - f) 
y oy T ' 


in the absence of external forces. This equation possesses a form similar to 
Eq. (4) . However, the latter eq'uation is more complex to treat since it 
includes teims involving Bf/Sc . Thus, one of the important extensions of the 
discrete ordinate method as applied to the present problem has been the treat- 
ment of derivatives in velocity space. Interestingly enough, the presence of 
external force terms in the Boltzmann equation would Introduce derivatives 
of this tsrpe so that knowledge gained in the present numerical solution for 
turbulence can be transferred back to rarefied gasdynamics. 

Since the moments required to compute flow properties of interest are 
obtained as Integrals over velocity space (cf. Eqs. (8)), it is the velocity 
variable which is chosen to be discretized. The set of points selected is 
denoted by {c^} j and a continuous function, say, g(y,Cy) is replaced by a 
set of functions gj^(y), o’ = 1,2, ... , S. The same procedure is applied to 
each of the dependent variables. The integrations over c to form moments may 
then be accomplished by numerical quadrature employing appropriate weighting 
functions. 
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CO 


(15) 


_S 

J 0(Cy) g(y, Cy) dCy s. ^ 0(c^) g^(y) 

-00 a=l 

where is a function of c , and are the weights in the quadrature. 

The details of the resulting equations are shown in Appendix A. 


Finite Difference. Methods 

Since derivatives with respect to both y and c^ are first order, the 

first approach taken in forming finite difference equations from the differential 

equations was to use simple forward or backward differences, depending upon the 

direction of integration. However, this first order scheme resulted in some 

numerical error in the region near the wall. This is illustrated in Figure 2 

for the solution obtained for Reynolds stress with the Sf/Sy = 0 boundary 

condition. It is ejected that c c should remain constant in the turbulent 

X y 

zone between the plates. As seen in this figure the deviation from a constant 
value is approximately one per cent. Although in most cases this is quite 
adequate in terms of accuracy, it was the non-constancy of the Reynolds Stress - 
as opposed to the absolute error - that was of some concern. The results for 
the Chapman-Enskog boundary condition show a similar, but exaggerated, behavior 
in that c^c^ varied approximately fifteen percent across the turbulent zone of 
the channel. Some of this variation is likely due to the model assumed for the 
boundary condition; however, it was deemed important to reduce numerical errors 
so that effects of physical modeling and numerical modeling could be more 
clearly delineated. 

Therefore, the possibility of employing a more accurate finite difference 
form was investigated. If a function f (y, c^) is expanded in a Taylor series 
about a point (y^, c^) , one may write 


^i-l,o 


= f. 







(Ay)^ - T + ... 


and 
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f . 


i-2,a 


= f. „ - 2 


i,a oy x.cr oy-^ x,a 


for a constant spacing Ay. Eliminating the second derivative terms, there 
results 


4f. T „ 




Solving for the first derivative gives 


(m 

\.3y/ 


3f. „ - 




f. 


i-l,g i-2,g 


ijO' 


2Ay 


+ 0 


[(Ay)^] 


Thus, this backward difference scheme has a truncation error of order (Ay) as 
compared to that of order (Ay) for the simple forward difference. A similar 
form for a second order forward difference scheme can be developed, yielding 


^ x,<J 2Ay 


(I). 


In the implicit finite difference scheme employed, the distribution 

function and its derivative with respect to y are to be evaluated at the same 

grid point. Consequently, the finite difference formulae shown above are more 

appropriate than the usual central difference scheme. 

A similar approach can be enployed for deriving expressions for 

fsf/Sc 1 . However, the spacing of discrete velocity points is necessarily 

^ i,cr 

variable so that efficient use of quadrature can be achieved. Therefore, it is 
preferable to obtain a second order finite difference expression from a 
Lagrange interpolation formula^^. This results in the following differences: 
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(backward) 


( forward) 


(if). 




y i,a ^°a-2 " i,0'-2 


, - ^a- 2 ^ 

^^a-1 “ %-2^^^cy-l " 


- ^g-l - °g-2 

(c - e _)(c -c ,) 
'' a g-2 g cr-l' 


+ 0 


K^O] 


(If). 


^^*^g “ *^g+l " °g+2^ 




y (Cjj - - c^+g) ’ 


^^g ~ ^g+ 2 ^ ^ 


^g+ 1 ^ 


^"^g+2 *^cr^^'^gH:2 '^cr+1^ 


^i,C7+2 



( 15 a) 


( 15 b) 
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As a consequence of the two-stream nature of the distribution functions, 
the choice of the difference scheme (either forward or backward) is readily 
prescribed. For the "positive" stream (c^ > O), the computations should proceed 
from + “ (where boundary conditions with respect to velocity space are known) 
to zero and from the lower boundaiy point (where conditions with respect, to 
physical space are known) to the upper boundary. Thus, the forward difference 
in velocity space and the backward difference in physical space are employed. 

For the "negative" stream (c^ O) the reverse is true. There, the integration 
proceeds from - “ to 0 in c and from upper boundary to lower boundary in y. 
Thus, the backward difference in velocity space and the forward difference in 
physical space are utilized. When these forms are substituted for the deriva- 
tive terms in Eqs. (9), a set of difference equations for the reduced distri- 
bution functions is obtained. These equations are given in detail in 
Appendix B. 

Iterative Scheme 

The resulting equations must be solved by an iteration process since they 
contain terms which depend upon the macroscopic properties. Therefore, initial 
guesses are made for u, U, and e (an initial profile for Reynolds stress is 
not required) . The equations for the "positive stream" are then solved from 
the boundary point up to the centerline, symmetry conditions are applied, and 
the "negative stream" is then computed from centerline to boundary point. 

This completes one iteration and yields the approximation to the reduced distri- 
bution functions. From these, new profiles for the macroscopic quantities are 
found and stored for use in the second Iteration. If integral constraints are 
required at the boimdary point, these are imposed before the second iteration is 
begun. These will be discussed in the next section. 

This iterative process continues until satisfactory convergence is obtained 
for the macroscopic properties. 

Constraints at Boundai-y 

The form of the distribution at the boundary point dictates the constraints 
which must be applied. 
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Zero Gradient Boundary Condition . In «i5)loyir!g the ho-uiKiary condition found 
from setting 9f/9y = 0, it is neeessaiy to specify U and (hence, the 
value of the wall shear stress) . This is the value of u^ that would result 
if the mean velocity were logarithmic between the plates. Equations (U) are 
then solved, subject to these constraints, to give the boundary conditions on 
g, j, h. These conditions are then fixed for all iterations. The values of 
u and du/dy at the boundaiy point are determined from the law of the wall, 
utilizing the assumed value of u^- 

Chapman -Enskog Boundary Conditions. As discussed ,the Chapman-Enskog form of 
the distribution function may be used as a boundary condition on the positive 
stream. With this form it is possible to deduce the value of wall shear from 
the solution to the equations, rather than requiring an a priori assumption on 
u^. This is achieved by applying appropriate integral constraints upon the 
outgoing and incoming streams at the boundary point. If the Chapman-Enskog 
forms are written for the reduced distribution functions, there results 


.( 1 ) = 


= G “ c [ -% (H + c^ G) - 5 g1 

Re^ dy y L y2 ^ y ^ J 


(l6a) 


,(i) = . 


Re* Ij2 dy 


(l6b) 


= H - ^ 






5H 


(l6c) 




+ U g 


( 1 ) 


(l6d) 
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Using these equations as houndary conditions on the outgoing or positive stream 
at the lower houndary point, the first iteration may hegin once initial guesses 
are posed for the macroscopic quantities. Then, upon marching hack from the 
centerline of symmetiy, certain quantities must he re-evaluated before the 
second Iteration can proceed. These are 


U 


V 

T ^ 
Ee^u3 dy 


and 


V 

T du _ 

Re^ dy X y 


In the present study the following constraints have heen applied: 


c c 
X y 


V 

- 1 - ^ =-10 
Re* dy 


(17a) 


J g dCy = Jg-dCy + J g- dCy = 1.0 


(17b) 


y 


= I 


' wmCO 


J 


c g do = 0 

y y 


(17c) 


The first of these states that the Reynolds stress at the boundary point is 
equal to the wall shear stress (the viscous stress could he included in this 
equation and, in fact, several calculations have heen performed over the course 
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of the study in which this has been done) . The second condition states that the 
probability of finding a fluid element with velocity between and is 

unity, while the third requires that the time average of the fluctuating 
vertical velocity component be zero. 

If the Chapman-Enskog forms of Eqs. (l6) are substituted into Eqs. (l 7 b,c), 
there results 


and 


V 

T 


Re^u3 


<3y 



(0.5 - c^) 
u 



c 


2 


where 


and 


°1 = 


- J 


g dc 


y 


"2 = 


c g dc 

y y 


(18a) 


(18b) 


The latter two integrals are computed numerically at the end of each iteration, 
based upon the current iterate for the g (or incoming stream) distribution. 
Thus, parameters in the outgoing stream may be readjusted at each iteration to 
conform with the imposed constraints. The value for u^ , and hence wall 
shear stress, is obtained by requiring that the computed value of u at the 
boundary point fit the logarithmic relation for law of the wall. 


u 


1 


K 



+ B 
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It is emphasized that this is the only point, under this Chapman-Enskog scheme 
for boundary conditions, at which law of the wall is assumed to hold; and this 
is implemented only to avoid using the statistical model for turbulence within 
the region where viscous stresses are comparable to or larger than Reynolds 
stresses. 



RESULTS 


All calculations reported here for the present study were made at a Reynolds 
number (Re = u^d/v) of 17,000. 

Zero-Gradient Boundary Condition 

The motivation for deriving this boundary condition and applying it to the 
Couette flow problem was twofold. First, it was important to determine whether, 
under appropriate assunptions, the statistical model equation could reproduce a 
turbulent flow within which the mean velocity varied logarithmically and the 
Reynolds stress remained constant. Since it is known from experiments that such 
a region exists near the wall for many turbulent flows, the capability of the 
model to recover this result is a logical first test of its validity. Second, 
it was expected that such a boundary condition might potentially be applied to 
more general situations than the Couette flow due to, as mentioned above, the 
existence of a limited logarithmic region for many boiuadary layer flows . 

So that the expression for e jji the statistical model equation is con- 
sistent with the law of the wall in the logarithmic region, the production is 

17 

set equal to the dissipation, giving 

e = 0.7968 U^/y (19) 

The zero gradient boundary conditions obtained from the numerical solution 
to Eqs. (11) and the equation for e , Eq. (19), were applied to the equations 
for Couette flow, Eqs. (9). The numerical scheme was that of second order finite 
differences derived in the previous section. An initial guess for mean velocity 
which corresponded to the law of the wall variation was assumed, and the iterative 
process was begun. Constants in the logarithmic mean velocity profile were taken 
as H = 0.4l and B = 5*0 . After 45 iterations, the mean velocity profile had 

converged to and remained within O.O5 percent of the logarithmic profile; and the 
dimensionless Reynolds stress was constant to three significant digits at - 0.992 
(exact value is - l.O) . The numerical procedure has been showa to give a unique 
solution for different' initial guesses for the velocity profile. For example, 
in one case the logarithmic profile was used as an initial guess and in another 
a linear profile for mean velocity was employed. The converged results agreed 
for both examples. 
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Solutions obtained usii^ the zero gradient boundary condition have clearly 
demonstrated that the Lundgren’s model equation is a reasonable one and have 
given encouragement that accurate results may be obtained from such a statistical 
approach. Further, it is believed that these results have demonstrated the 
numerical accuracy of the discrete ordinate and difference schemes presently 
employed. 

Chapman-Enskog Boundary Condition 

Although it is established that a logarithmic region exists near the wall 

for many boundary layer flows, this does not imply that such a region will 

extend across the entire field for the Couette flow case. In fact, experimental 

12 18 

data indicate that such may not be the situation ’ . Although Johnson's 

data for Couette flow can be made to fit the law of the wall if h = 0.4115 
and B = 5*6, the fact that these "constants" may not fit other experimental data 

13 

indicates a lack of rigor in specifying a unique resxilt. Furthermore, Eeichardt's 
data for mean velocity do not fit a logarithmic variation very well for any pair 
of these constants. 

If the experimental value of u^ quoted by Reichardt is employed and 
reasonable values of h and B are assumed, the logarithmic velocity profile 
does not pass through zero at the centerline and does not fit the experimental 
data well. On the other hand, by selecting a rather large value of B (7.456) 
the logarithmic profile can be made to satisfy zero velocity at the centerline, but 
this does not follow the experimental data elsewhere. 

The Chapman-Enskog form of the distribution function for the outgoing 
stream was employed as a boundary condition in a series of calculations. The 
first results reported here were obtained with the first order difference scheme. 
Calculations utilizing the improved second order scheme are currently in pro- 
gress and some results will be given. The major difference obtained by using 
these two numerical methods is improved accuracy near the boundary with the 
second order technique. Numerical solutions have been calculated using both 
the algebraic (Eq. (9e) ) and differential (Eq. (9f) ) equations for the dissipa- 
tion e . In all cases reported, the solutions were unique and convergent. 

This was tested by employing different initial profiles, various grid spacings, 
and different boundary point locations. 
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Figure 3 illtistrates the results for mean velocity. The ejqperimental 
data of Relchardt are plotted for comparison. The logarithmic profile corre- 
sponding to K = 0.4l and B = 5.0 is also shown. Both of the solutions 
using the statistical model equation fall very near to the eicperimental data. 

The case for e ^ u^/y is particularly close. Some differences can be observed 
between the solution obtained assuming that dissipation equals production and 
that found when the differential equation for g is utilized. These differences 
are relatively minor. Chung's solutions are also plotted for comparative 
purposes. The ratio of friction velocity to plate velocity quoted by Reichardt 
was 


— ^ 0.0425 (Re = 17,000) 
w 


while the calculated values were 


and 


— = 0.04487 


(Using Eq. (9e) for e ; Re = 17,000) 


% 

— = 0.04467 (Using Eq. (9f) for s j Re = 17,000) 
w 


Chung's calculated values^ were 0.03284 (Re = 365OOO) and 0.0405 (Re - 9j6lO). 

Resixlts for Reynolds stress are given in Figure 4. The ejcperimental 
1 '^ 1 

data'^'^’ are shown for comparison. The calculated results show a variable 
Reynolds stress which decreases near the wall boundary. The decrease is not a 
consequence of increasing viscous stress but rather is a result of some numerical 
error with the first order scheme and of modelling the boundary condition tising 
the Chapman-Enskog form. Despite this variation, the calculated results are 
within 1 15 per cent of the experimental values. Experience with the second 
order difference scheme shows substantial Improvement in the Reynolds stress 
profile near the wall. Again, Chung's results^ are shown for comparison. 
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The profiles for kinetic energy are given on Figure 5* There is very 
little difference between the results for the two models for e , and both 
profiles show a slight decrease in turbulence kinetic energy as the wall 
boundary point is approached. 

An interesting aspect of the statistical approach to turbulence is that 
the distribution functions may be calculated. This may offer insight into 
turbulence mechanisms and aid turbulence modeling. Figure 6 illustrates the 
g distribution employed as a boundary condition for both the zero gradient and 
Chapman-Enskog cases, while Figure 7 shows a similar graph for j. (it is 
emphasized that only the outgoing st3ream is modeled with these functions). The 
Chapman-Enskog distributions are somiewhat broader than those for 3f/3y = 0. 

This is particularly noticeable for the j distribution function. It should 
be pointed out again that a Gaussian distribution would lead to j = 0 for all 
c^ and to a zero value of Reynolds stress. It should also be pointed out 
here that the Chapmian-Enskog distribution does not satisfy the governing differ- 
ential eq-uation at the boundary. Therefore some gradients can be expected with 
such an approach. 

Figure 8 shows results for g at several y locations across the channel 
for the Chapman-Enskog boundary condition case. The two-stream nature of the 
distribution function when employing this boxmdary condition is clearly evident. 
Near the wall (y = 0.125) there is a noticeable discontinuity in g at c =0. 
This discontinuity gradually decreases as y Increases. This is a consequence 
of the influence of the relaxation term in the governing equation. Physically, 
the fluid elements are interacting to smooth out the distribution function. 

Recently obtained results employing the Chapman-Enskog boundary condition, 
the differential equation model for e , and the second order difference scheme 
are shown in Figures 9 throx:igh 13 . Figure 9 illustrates the results for the mean 
velocity profile and compares these with Reichardt’s experimental data. The 
comparison is quite favorable. The Reynolds stress is given in Figure 10. It 
can be seen that, although the calculations do not give a flat profile, the sharp 


decline in c c near the wall experienced with the first order scheme is sub- 
X y 

stantially reduced by employing the second order difference form. The calculated 
value for is 0.0^^4377, compared with the value of 0.0ij25 deduced from 

Reichardt’s data. The turbulence kinetic energy, shown in Figure 11, is also 
reasonably constant across the channel. Figure 12 Illustrates the results obtained 


29 


I 


for e from Equations (9e) and (9f)» The two solutions compare quite well. 

This indicates that the algebraic expression for e , Eq. (9e) , is very good for 
Couette flow. 

Another interesting feature of the present approach is that it is also 

possible to compute the contribution of the y-component of velocity fluctuations 

to the turbulence kinetic energy. This contribution, as shown in Figure 13, is 

almost constant except for a slight variation near the boundary point. In 

T 

regions away from the boundary point, the solution obtained shows that c is 
about 89 per cent of U^. In this region, as seen from Figure 8, the distribu- 
tion function, g , does not vary with Y. For such cases, it can be shown from 
Equation (9a) that 


2 

c 

y 



For the value of K used in this region, this ratio is about O. 89 . This ratio, 
however, is quite different from the isotropic value. Thios, the results obtained 
with the inproved difference scheme are quite encouraging. 
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coucaajsiONS 


The present study has accomplished the following: 

(i) A convergent numerical scheme employing a 

combination of the discrete ordinate method and 
finite differences has been developed for solving 
the one -dimensional form of Lundgren's model 
equation for turbulence. 

(ii) Physically realistic boundary conditions for the 

distribution function and models for the tiurbulence 
dissipation rate have been examined. 

(iii) Lundgren’s equation has been proven to yield rea- 
sonable results for mean velocity, Reynolds stress, 
and turb'ulence kinetic energy for the case of simple 
Couette flow. 

Thus, it is believed that this research has yielded important contributions to 
the understanding and modelling of tTrrbulent flows; and further, that the 
knowledge gained provides a basis for additional studies in fundamental aspects 
of turbulence. 

The comparisons of theory and experiment reported here indicate that the 
statistical approach taken by L-undgren provides an accurate description of a 
simple case of wall-bounded turbulence — Couette flow with no pressure gradient. 
However, there are only limited experimental data available for comparison. An 
appropriate extension of the present work would be to consider the case of 
channel flow, for which more extensive measurements are published. One of the 
primary areas of study should be the use of such data to better model the 
boundary conditions for the distribution function. Further refinements in the 
statistical model itself should be considered if comparisons of theoretical and 
experimental results indicate difficulties with the basic BGK-type approach 
used by Lundgren. Based on the experience with Couette flow, numerical solu- 
tions should be relatively straightforward and inexpensive (in terms of computer 
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time) for one -dimensional problems -which include pressure gradient and 
chemical reactions. 

Computation time to achieve a converged solution for the zero-gradient 
bo-undary condition was about 12 minutes on a UKIVAC U08. For the Chapman- 
Enskog boundary condition, total computation time varied from 25 to 50 minutes, 
depending on the choice of initial profiles. 

Althoiigh the n-umerical techniques would be more complex and time-consuming, 
certain simple two-dimensional problems could also be attacked using the model 
equation and solution method employed in the present research. Examples of 
this are free shear layers and boundary layer flows. However, it is thought 
that the real value in employing the statistical approach to t-urbulent flows 
examined here is in furthering basic understanding, as opposed to developing 
a practical computational tool. 
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APPEMDIX A 


The governing differential equations, Eqs. ( 9 ), contain terms which are 
defined as moments of the distribution functions. Since these moments are 
obtained as integrals over velocity space, as shown in Eqs. (8), it is 
convenient to discretize the velocity variable. The set of discrete velocity 
points selected is denoted by {c^} , and a continuous function, say 
g(y, Cy), is replaced by a set of functions gg.(y), CJ = 1, 2 . . . , S. This 
procedure is applied for each of the dependent variables. To evaluate the 
required moments, numerical quadratures employing appropriate weighting 
functions may be used. These are of the form 



g(yj Cy) 


S 

dCy 2 0(c^) g^(y) % 

a=l 


(A-1) 


where ^ is a function of c , and W are the weights in the quadrature. 

y 

The specific values for the discrete ordinates, c^ , depend upon the 
quadrature selected. In the present study an eleven-point closed-type 
Eewton-Cotes quadrature is employed for Integration. This quadrature requires 
equally spaced points in the interval of integration. To reduce the computing 
time, the interval of Integration is divided into many sub -intervals and 
appropriate spacings are chosen in each sub-interval. 

After some numerical experimentation, a total of discrete velocity 
points were utilized to insure sufficient accuracy in the integrations. 

When discrete ordinates are employed, the governing equation {3b-), then 
becomes a system of differential equations of the form 




<3y 




3u2 Of) } 

y cr 


(A-2) 
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•where a = 1, 2, . . • , S. In this system, is a ftuiction of y only, 

and it represents the function g(y, c ) e’va.l’uated at the discrete velocity 

y 


point Cy = c^. 
(A-2), the term- 


In the n'umerical procedure used to solve the set of equations 



is replaced hy a finite difference approximation as 


described in Appendix B. Thus, the governing partial differential equation, 

Eq, (9a) is approximated by a set of ordinary differential eq'uations of the 
form sho'wn in Eq. (A-2) . A similar set of equations is obtained for each of 
the reduced distribution functions. These sets are then solved for the various 


discretized functions, and the results are used in the numerical quadrat’ures 
representing integrations over velocity space (Eq. (A-l)). 
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APPEHDIX B 


The non-dimens ionalized governing equations for the reduced distribution 
functions are 


and 


y Sy 




(B-1) 


■ 7 ri - - “y f) S 


(B-2) 


+ _1_ c ^ 

3if 


Bh 

V 9y 


i (h - h) 

+ 2 (:^ 

,2 

^ - P 

- c 

T \ / 

\Re* 

dv^ ^ 

y <ay / 


(B-3) 
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(B-4) 
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Using the second-order finite difference schemes outlined in the 
section UUMEEIGAL APPROACH, the Equations (B-1) - (b- 4) can he approximated 
for each node (i,cr) as shorn below. The subscript "1" denotes the 1*^'^ point 
in the physical space, y^ , and the subscript "ct" represents the discrete 
velocity point c^^. 

finite -difference equations for "positive stream", (c^ > O) can be 
obtained by using backward differences in physical space and forward differences 
in velocity space. The reduced distribution functions for "positive stream" 
are denoted by a superscript "+". Thus, Equation (B-1) can be written, for 
> 0, as 


where 


2(^y) ~ ®i-2,aj ~ ? l^i,a " H,o) 


^ i + ^ i ^CT / ++ ++ + + \ 

^ ®i,a 2 W ®i,CT ^2 ®i,a+l ®i,o+2/ 


3U, 


+ ^ ~ %+l “ ^^ 0 + 2 ^ 

Ut 


" ( 




= (c. 


K - ^0+2^ 


0+1 “ ^'^ 0+1 “ ‘^ 0 + 2 ^ 
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and 






^^a+2 ~ ~ ‘^CT+l^ 


Solving for g^ , one gets 


g+ = { ^u£ + ffS.') (2 g+ - - ) 

^i,a I T \AyJ \ ^i-l,a 2 ®i-2,a/ 


(B-5) 


® -c_ 


, 'i O' ( + + , T,'^ + 

+ ^ \^2 Si,o+i 

3U. 


,a+l ^3 ^i,a+ 2 . 


:)} 





Similar expressions for positive stream can be derived for the other reduced 
distribution functions. 

Finite-difference equations for "negative stream", (c^ < O) could be 
derived by using fonrard differences in physical and backward differences in 
velocity space. Equation (B-1) can then be written for "negati-ve stream", 
c^ < 0 ) indicated by a superscript "-", as 
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2(Ay) 


®i+l,a " ®i+2,CT " 


3U, 


2 ®i,a 2 (°1 ®l,a-2 '*' °2 ®i,a-l '*' °3 ®i,c) 


where 


\ = 


(c - c ) 
a a-1 


^^a-2 " ‘^a-1^ *‘‘^a-2 


and 


- =a-2) 


■ '^a-2^^°cj-i " 


(2c - c , - c „) 

g CT-1 g-2 

(^g " V2^(^g " W) 


Solving for g. „ > one gets 


= 1 _ia£ + — Sl_ (i g“ - 2 g“ ) 


(B-6) 


^g (°1 %,g-2 '*' °2 ®i,g-l 
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Similar easpressions for negative stream can te derived for the other reduced 
distribution functions. 
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